The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48_1/routh

Internet Message Format | 1995-03-31 | 13KB

From helens!shelby!rutgers!uwm.edu!ux1.cso.uiuc.edu!iuvax!noose.ecn.purdue.edu!en.ecn.purdue.edu!ryoder Mon Jul  2 12:08:13 PDT 1990
Status: RO

Article 2050 of comp.sys.handhelds:
Path: helens!shelby!rutgers!uwm.edu!ux1.cso.uiuc.edu!iuvax!noose.ecn.purdue.edu!en.ecn.purdue.edu!ryoder
>From: [e-mail] (Robert W Yoder)
Newsgroups: comp.sys.handhelds
Subject: Routh-Hurwitz
Keywords: Routh-Hurwitz
Message-ID: <[e-mail]>
Date: 2 Jul 90 16:36:52 GMT
Organization: Purdue University Engineering Computer Network
Lines: 245

The following is the response I received from Paul Dujmich after I sent him a
Routh-Hurwitz program I had written.  I had never posted it because I had not
gotten around to adding the necessary code to cover the special cases.  He
expanded it to cover those cases, so here it is.

Robert Yoder
[e-mail]

--------------------------------------------------------------------------------
Bob

I took the liberty of modifying your basic Routh program for the purpose of
adding special case checking. Most of the special case code is in a separate
subroutine I call "ARYCHK". It was just easier to do it this way.
The new code checks for left-most column zero, as well as zeroes all the
way across a row (premature termination). I believe these two are the
only special cases to worry about. I have tested the new version on my
homework problems, which I had worked out by hand. They included all the
special case conditions mentioned above. The program produced accurate
results for every polynomial that was tried.
One note is in order.
For left-most zero, we use the epsilon method, replacing the 0 with epsilon,
a very small number approaching zero.Then we take the limit of the term,
as epsilon approaches 0. In the program, I use  +/- 1E-50 to represent
epsilon. This works fine, but produces slightly different numbers in the
finished array (compared to the limit method). But, we really don't care
about the actual values anyway, all we care about is the sign changes, and
the program has been correct in every case tested. For zeroes all the way
across a row, I do the accepted method of forming a new Aux polynomial
from the coefficients in the row just above the row of zeroes. Then I take
the derivative of the Aux polynomial, and use the new coefficients to
replace the row of zeroes. This has also worked flawlessly for every thing
that was tried.

Please feel free to post this version on net.handhelds if you want. I'm sure
there are other EE/BSET people who would like to have it. Thanks again for
your reply to my posting, and for your original program. It sure was a lot
easier not having to start from scratch. Well anyway, here's the updated
version.


Thanks Again

Paul
[e-mail]

============================================================================
%%HP: T(3)A(R)F(.);
@ PROGRAM NAME:  ROUTH
@ DATE:          6-30-90
@ PURPOSE:       Does Routh-Hurwitz Stability Testing
@ ARGUMENTS:     1. An array in level 2 containing coefficients of the
@                   characteristic polynomial from the transfer function.
@                2. The ORDER of the characteristic polynomial in level 1
@ NOTE: This program requires "ARYCHK", a subroutine that checks for special
@       cases, to operate.


@Routh-Hurwitz  Program
@This program tests the characteristic polynomial of a transfer function
@for stability. The program expects an array as input, which is made up of
@the coefficients of the characteristic polynomial. The array must be in stack
@level 2 when the program is executed. The program also expects the ORDER of
@the characteristic polynomial in level 1.
@The program outputs a Routh-Hurwitz array as it's output.
@The stability of the system being tested can be obtained
@by examining the left-most column of the output array. If all elements in the
@left-most column of the array are of the same sign, the system is STABLE, and
@all of the polynomial roots are in the left-half plane of the complex plane.
@If there are any sign changes within these elements, then the system is
@UNSTABLE, and the polynomial has roots in the RIGHT HALF PLANE of the complex
@plane, or
@it has roots on the IMAGINARY AXIS of the complex plane. The NUMBER of sign
@changes within the left-most array elements give the number of roots in the
@right half plane. The program tests each input array for the following
@SPECIAL CASES, and applies a fix, so that the array may be completed:
@            Special Case 1. Zero in left-most array element
@            Special Case 2. Zeroes all the way across a row


\-> ord             @load in polynomial ORDER

\<< DUP             @make copy of input array on stack
    SIZE            @return dimensions of input array {2 2} on stack
    LIST\->         @dump the array dimensions from list to stack
                    @ 2   (rows)
                    @ 2   (columns)
                    @ 2   (there were 2 elements in the list)
    DROP            @throw away the #elements in list, don't need it
    'cols' STO      @save # columns of input array
    DROP            @drop off # columns from stack
    cols 3 *        @new # of rows = 3* original # columns
    cols            @orig number of columns
    2 \->LIST       @put last 2 elements in a list
    RDM             @redimension input array to size specified by list
    'a' STO         @save the redimensioned array
    3 'r' STO       @start calculating at row 3
    1 CF            @clear flag 1 in case it is set

@**************** calculate # of coefficients in orig array  *************

ord 1 +             @increment polynomial order to get # coefficients
'coef' STO          @save it

@************** Check row 2 of input array for special cases ***************

CLEAR                 @clear stack
3 'r' STO
ARYCHK                @check row 2 of the input array for special cases
                      @if they exist-fix them.

@********************** main DO LOOP ****************************************
3 'r' STO           @start calculating at row 3
  DO
     1 cols 1 -     @set up FOR loop for 1 less than number of columns
     FOR c
         a          @put array on stack
         r          @put current row on stack
         c          @put current column within current row on stack
         2 \->LIST  @make a list of last 2

         @************* array element calculation ************************
         '(a(r-1,1)*a(r-2,c+1)-a(r-2,1)*a(r-1,c+1))/a(r-1,1)'
         EVAL
         PUT        @put the calculated element back into the array
         'a' STO    @save the updated array
     NEXT
     1 'r' STO+     @increment to next row
     ARYCHK         @check for special cases
UNTIL 1 FS?         @run DO loop until flag 1 is set
END
a                   @ put array on stack
r 1 - cols 2 \->LIST
RDM                 @ redimension the array
{ a cols r s cnt coef hc pwr}  @list of variables to nuke
PURGE               @nuke them
1 CF                @clear flag 1
\>>
=============================================================================
%%HP: T(3)A(R)F(.);
@ PROGRAM NAME:  ARYCHK  (subroutine for Routh-Hurwitz Program)
@ DATE:          6-30-90
@ PURPOSE:       checks a row of a Routh-Hurwitz array for special cases
@ARGUMENTS:      Variable 'r' must exist and contain the row of the Routh-
                 @Hurwitz array (+1) to be checked.

\<<
     IF 'a(r-1,1)==0'                     @if leftmost element of last row == 0
               THEN
                   IF 'r-1 > coef'        @normal valid array termination
                           THEN 1 SF          @set the flag

                          @special case checking
                            ELSE 0 'cnt' STO        @clear cnt
                            1 cols                  @FOR loop goes from 1 to
                                                    @ # columns
                            FOR c
                                'a(r-1,c)' EVAL     @count # valid array
                                                    @elements in row where 0
                                                    @occurred
                                0 \=/               @element !=0?
                                'cnt' STO+          @if yes, then valid
                            NEXT
                            IF 'cnt==0'             @the whole row had 0's

                               @premature termination of array (all 0's in row)
                                  THEN
                                     coef r 2 -  -  @calculate first power of
                                                    @S for new polynomial
                                     'pwr' STO      @save it
                                     1 cols         @FOR loop goes from 1 to
                                                    @width of a column
                                     FOR c
                                         'a(r-2,c)'
                                         EVAL       @get term coefficient
                                         "*S^"
                                         +          @add it to string
                                         pwr        @get the exponent
                                         +          @add it to string
                                         "'"
                                         SWAP
                                         +          @make algebraic object
                                                    @ inside a string
                                         OBJ\->     @remove string quotes
                                                    @object is now algebraic
                                         1 'S' STO
                                         'S'        @variable of differentiation
                                         \.d        @take the derivative
                                         EVAL       @evaluate at S = 1
                                         a SWAP     @put array on stack
                                         r 1 - c 2 \->LIST @elem to change
                                         SWAP
                                         PUT        @correct the element
                                         'a' STO    @save updated array
                                         pwr 2 -    @decrement exponent by 2
                                         'pwr' STO  @save new power
                                         IF pwr 0 < @if exponent < 0 then
                                                    @terminate the FOR loop
                                            THEN cols 1 +
                                                 'c' STO
                                         END
                                     NEXT
                                     'S' PURGE      @ done with S for now

                              @zero in left-most column only
                                ELSE 0 'cnt' STO    @clear cnt
                                     1 coef         @FOR loop goes from 1 to
                                                    @ # of coefficients
                                     FOR c
                                          'a(c,1)'   @count the number of left
                                          EVAL       @most elements that are
                                          0 <        @negative
                                          'cnt' STO+ @cnt has total
                                     NEXT
                                     coef 2 /        @find half #coefficients
                                     IP              @only the integer part
                                     'hc' STO        @save it in hc
                                     IF 'cnt > hc'   @more than 1/2 the
                                                     @coefficients are negative
                                       THEN  a       @get array
                                             r 1 - 1 2
                                             \->LIST    @array element to change
                                             -1.E-50    @make it neg
                                             PUT        @put in array
                                             'a' STO    @save changed array
                                       ELSE  a          @get array
                                             r 1 - 1 2  @select element
                                             \->LIST
                                             1.E-50   @make it pos
                                             PUT      @put in array
                                             'a' STO  @save changed array
                                     END

                            END
                   END

     END
   \>>
==============================================================================

--
--------------------------------------------------------------------------------
Robert Yoder

Internet: [e-mail]